packages <- c('sf', 'tidyverse', 'tmap')
for (p in packages) {
if (!require(p, character.only = T)) {
install.packages(p)
}
library(p, character.only = T)
}
## Loading required package: sf
## Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
## Loading required package: tidyverse
## ── Attaching packages ─────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.3 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 2.0.0 ✓ forcats 0.5.1
## ── Conflicts ────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## Loading required package: tmap
setwd("~/IS415_blog/_posts/2021-09-04-take-home-exercise-1")
library(readxl)
covid_monthly <- read_excel("data/aspatial/covid_monthly.xlsx")
jarkarta <- st_read("./data/geospatial/BATAS DESA DESEMBER 2019 DUKCAPIL DKI JAKARTA/BATAS_DESA_DESEMBER_2019_DUKCAPIL_DKI_JAKARTA.shp")
## Reading layer `BATAS_DESA_DESEMBER_2019_DUKCAPIL_DKI_JAKARTA' from data source `/Users/pengtaixu/IS415_blog/_posts/2021-09-04-take-home-exercise-1/data/geospatial/BATAS DESA DESEMBER 2019 DUKCAPIL DKI JAKARTA/BATAS_DESA_DESEMBER_2019_DUKCAPIL_DKI_JAKARTA.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 269 features and 161 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 106.3831 ymin: -6.370815 xmax: 106.9728 ymax: -5.184322
## Geodetic CRS: WGS 84
jarkarta_sel = jarkarta[,1:9]
jarkarta_main <- jarkarta_sel[!grepl("PULAU", jarkarta_sel$DESA),]
plot(st_geometry(jarkarta_main))
names(covid_monthly)[names(covid_monthly) == 'ID_KEL'] <- 'KODE_DESA'
jarkarta_covid <- merge(x = covid_monthly, y = jarkarta_main, by = "KODE_DESA", all.x = TRUE)
jarkarta_covid <- st_as_sf(jarkarta_covid)
jarkarta_covid <- jarkarta_covid[!is.na(jarkarta_covid$OBJECT_ID),]
jarkarta_covid <- jarkarta_covid %>%
mutate(CASE_RATE = POSITIF / (JUMLAH_PEN / 10000))
cumulative_death_by_month <- aggregate(jarkarta_covid$CASE_RATE, by=list(Category=jarkarta_covid$Date), FUN=mean)
names(cumulative_death_by_month) <- c("Date", "Cumulative Case Rates")
print(cumulative_death_by_month)
## Date Cumulative Case Rates
## 1 2020-03-31 0.7608927
## 2 2020-04-30 3.3702865
## 3 2020-05-31 5.4575059
## 4 2020-06-30 8.7275655
## 5 2020-07-31 14.0024393
## 6 2020-08-31 25.2031728
## 7 2020-09-30 51.5304429
## 8 2020-10-31 76.9717084
## 9 2020-11-30 103.0931858
## 10 2020-12-31 142.6647998
## 11 2021-01-31 212.3110015
## 12 2021-02-27 278.2796102
## 13 2021-03-31 318.4017212
## 14 2021-04-30 318.4017212
## 15 2021-05-31 362.2194383
## 16 2021-06-30 465.2808297
## 17 2021-07-31 705.3066398
## 18 2021-08-31 705.3066398
jarkarta_covid <- jarkarta_covid %>%
mutate(DEATH_RATE = Meninggal / (JUMLAH_PEN / 10000))
cumulative_death_by_month <- aggregate(jarkarta_covid$DEATH_RATE, by=list(Category=jarkarta_covid$Date), FUN=mean)
names(cumulative_death_by_month) <- c("Date", "Cumulative Death Rates")
print(cumulative_death_by_month)
## Date Cumulative Death Rates
## 1 2020-03-31 0.08094349
## 2 2020-04-30 0.34819827
## 3 2020-05-31 0.46113478
## 4 2020-06-30 0.55789903
## 5 2020-07-31 0.69121319
## 6 2020-08-31 0.88272831
## 7 2020-09-30 1.21924179
## 8 2020-10-31 1.68455456
## 9 2020-11-30 2.02190139
## 10 2020-12-31 2.58311191
## 11 2021-01-31 3.50570280
## 12 2021-02-27 4.59829464
## 13 2021-03-31 5.43802194
## 14 2021-04-30 5.43802194
## 15 2021-05-31 6.35020052
## 16 2021-06-30 7.35146861
## 17 2021-07-31 10.60192369
## 18 2021-08-31 10.60192369
tm_shape(jarkarta_covid) +
tm_polygons("CASE_RATE") +
tm_facets(by = "Date")
tm_shape(jarkarta_covid) +
tm_polygons("DEATH_RATE") +
tm_facets(by = "Date")
get.var <- function(vname,df) {
v <- df[vname] %>% st_set_geometry(NULL)
v <- unname(v[,1])
return(v)
}
# percent map
percentmap <- function(vnam, df, legtitle=NA, mtitle="Percentile Map"){
percent <- c(0, .01, .1, .5, .9, .99, 1)
var <- get.var(vnam,df)
bperc <- quantile(var, percent)
tm_shape(df) +
tm_polygons() +
tm_shape(df) +
tm_fill(vnam,title=legtitle,
breaks=bperc,
palette="Blues",
labels = c("< 1%", "1% - 10%", "10% - 50%", "50% - 90%", "90% - 99%", "> 99%")) +
tm_borders() +
tm_layout(title = mtitle,
title.position = c("right", "bottom")) +
tm_facets(by = "Date")
}
percentmap("CASE_RATE", jarkarta_covid)
percentmap("DEATH_RATE", jarkarta_covid)
Distill is a publication format for scientific and technical writing, native to the web.
Learn more about using Distill at https://rstudio.github.io/distill.